home *** CD-ROM | disk | FTP | other *** search
/ Mac Magazin/MacEasy 12 / Mac Magazin and MacEasy Magazine CD - Issue 12.iso / Sharewarebibliothek / Anwendungen / Wissenschaft & Technik / Yorick / yorick11-ppc folder / testm.i < prev    next >
Text File  |  1995-06-28  |  6KB  |  217 lines
  1. /*
  2.    TESTM.I
  3.    Certify functions in fft.i and matrix.i
  4.  
  5.    $Id$
  6.  */
  7. /*    Copyright (c) 1994.  The Regents of the University of California.
  8.                     All rights reserved.  */
  9.  
  10. func testm
  11. {
  12.   elapsed= old_elapsed= [0., 0., 0.];
  13.   write, "testing fft routines...";
  14.   timer, old_elapsed;
  15.   fft_test, 95;
  16.   fft_test, 32;
  17.   fft_test, 12;
  18.   timer, elapsed;
  19.   timer_print, "FFT elapsed time", elapsed-old_elapsed;
  20.   write, "testing LU decomposition routine...";
  21.   timer, old_elapsed;
  22.   testLU, 75;
  23.   testLU, 16;
  24.   timer, elapsed;
  25.   timer_print, "LU elapsed time", elapsed-old_elapsed;
  26.   if( !is_void(_dgelss) && !is_void(_dgesvx) && !is_void(_dgelx) ) {
  27.     write, "testing QR decomposition routine...";
  28.     timer, old_elapsed;
  29.     testQR, 75, 75;
  30.     testQR, 16, 16;
  31.     testQR, 25, 21;
  32.     testQR, 21, 25;
  33.     timer, elapsed;
  34.     timer_print, "QR elapsed time", elapsed-old_elapsed;
  35.     write, "testing SVD routine...";
  36.     timer, old_elapsed;
  37.     testSVD, 75, 75;
  38.     testSVD, 16, 16;
  39.     testSVD, 25, 21;
  40.     testSVD, 21, 25;
  41.     testSV, 21;
  42.     testSV, 64;
  43.     timer, elapsed;
  44.     timer_print, "SVD elapsed time", elapsed-old_elapsed;
  45.   }
  46. }
  47.  
  48. func fft_test(n)
  49. {
  50.   index= 2*pi*(indgen(n)-1.0)/n;
  51.   z= sin(index*3);
  52.   zf= fft(z, 1);
  53.   z3= z2= array(0i, n);
  54.   z3(4)= -0.5i*n;
  55.   z3(-2)= 0.5i*n;
  56.   zb= fft(z, -1);
  57.   if (max(abs(zf-z3))>1.e-12*n || max(abs(zb-conj(z3)))>1.e-12*n)
  58.     write, "***WARNING*** failed 1D fft test";
  59.   if (n<=96 || force2d) {
  60.     z*= cos(index*2)(-,);
  61.     zf= fft(z, [0, 1]);
  62.     z2(3)= z2(-1)= 0.5*n;
  63.     zb= fft(z, [], [-1, 0]);
  64.     if (max(abs(zf-sin(index*3)*z2(-,)))>1.e-12*n ||
  65.     max(abs(zb-conj(z3)*cos(index*2)(-,)))>1.e-12*n)
  66.       write, "***WARNING*** failed first 2D fft test";
  67.     zf= fft(z, 1);
  68.     zb= fft(z, -1);
  69.     if (max(abs(zf-z3*z2(-,)))>1.e-12*n ||
  70.     max(abs(zb-conj(z3)*z2(-,)))>1.e-12*n)
  71.       write, "***WARNING*** failed second 2D fft test";
  72.   }
  73. }
  74.  
  75. func TDcheck(c, d, e, b, x, s)
  76. {
  77.   check= _(   d(1)*x(1)    +    e(1)*x(2),
  78.        c(1:-1)*x(1:-2) + d(2:-1)*x(2:-1) + e(2:0)*x(3:0),
  79.                             c(0)*x(-1)   +   d(0)*x(0)   );
  80.   if (max(abs(check-b))>1.e-9*max(abs(b))) {
  81.     write, "***WARNING*** "+s+" tridiagonal solution doesn't check";
  82.     write, "   max relative error is "+pr1((max(abs(check-b)))/max(abs(b)));
  83.   }
  84. }
  85.  
  86. func testTD(n)
  87. {
  88.   c= random(n-1);
  89.   d= random(n);
  90.   e= random(n-1);
  91.   b= random(n);
  92.   TDcheck,c,d,e,b,TDsolve(c,d,e,b), "1D";
  93.   b2= random(n);
  94.   x= TDsolve(c,d,e,[b,b2])
  95.   TDcheck,c,d,e,b, x(,1), "2D(1)";
  96.   TDcheck,c,d,e,b2, x(,2), "2D(2)";
  97.   x= TDsolve(c,d,e,transpose([b,b2]), which=2)
  98.   TDcheck,c,d,e,b, x(1,), "2D(1)/which";
  99.   TDcheck,c,d,e,b2, x(2,), "2D(2)/which";
  100. }
  101.  
  102. func LUcheck(a, b, x, s)
  103. {
  104.   check= a(,+)*x(+);
  105.   if (max(abs(check-b))>1.e-9*max(abs(b))) {
  106.     write, "***WARNING*** "+s+" LUsolve solution doesn't check";
  107.     write, "   max relative error is "+pr1((max(abs(check-b)))/max(abs(b)));
  108.   }
  109. }
  110.  
  111. func testLU(n)
  112. {
  113.   a= random(n,n);
  114.   b= random(n);
  115.   LUcheck,a,b,LUsolve(a,b), "1D";
  116.   b2= random(n);
  117.   x= LUsolve(a,[b,b2])
  118.   LUcheck,a,b, x(,1), "2D(1)";
  119.   LUcheck,a,b2, x(,2), "2D(2)";
  120.   x= LUsolve(transpose(a),[b,b2])
  121.   LUcheck,transpose(a),b, x(,1), "t2D(1)";
  122.   LUcheck,transpose(a),b2, x(,2), "t2D(2)";
  123.   x= LUsolve(a,transpose([b,b2]), which=2)
  124.   LUcheck,a,b, x(1,), "2D(1)/which";
  125.   LUcheck,a,b2, x(2,), "2D(2)/which";
  126.   x= LUsolve(transpose(a),transpose([b,b2]), which=2)
  127.   LUcheck,transpose(a),b, x(1,), "t2D(1)/which";
  128.   LUcheck,transpose(a),b2, x(2,), "t2D(2)/which";
  129.   ai= LUsolve(a);
  130.   err= max(abs(ai(,+)*a(+,)-unit(n)))/max(abs(ai));
  131.   if (err>1.e-9) {
  132.     write, "***WARNING*** LUsolve inverse is fishy";
  133.     write, "   max relative error is "+pr1(err);
  134.   }
  135. }
  136.  
  137. func QRcheck(a, b, x, s)
  138. {
  139.   check= a(,+)*x(+);
  140.   err= max(abs(check-b))/max(abs(b));
  141.   if (err>1.e-9 && dimsof(a)(2)<=dimsof(a)(3)) {
  142.     write, "***WARNING*** "+s+" QRsolve solution doesn't check";
  143.     write, "   max relative error is "+pr1(err);
  144.   }
  145. }
  146.  
  147. func testQR(m,n)
  148. {
  149.   a= random(m,n);
  150.   b= random(m);
  151.   QRcheck,a,b,QRsolve(a,b), "1D";
  152.   b2= random(m);
  153.   x= QRsolve(a,[b,b2])
  154.   QRcheck,a,b, x(,1), "2D(1)";
  155.   QRcheck,a,b2, x(,2), "2D(2)";
  156.   x= QRsolve(a,transpose([b,b2]), which=2)
  157.   QRcheck,a,b, x(1,), "2D(1)/which";
  158.   QRcheck,a,b2, x(2,), "2D(2)/which";
  159. }
  160.  
  161. func SVcheck(a, b, x, s)
  162. {
  163.   check= a(,+)*x(+);
  164.   err= max(abs(check-b))/max(abs(b));
  165.   if (err>1.e-9) {
  166.     write, "***WARNING*** "+s+" SVsolve solution doesn't check";
  167.     write, "   max relative error is "+pr1(err);
  168.   }
  169. }
  170.  
  171. func testSVD(m, n)
  172. {
  173.   a= random(m,n);
  174.   s= SVdec(a, u, v);
  175.   if (anyof(s(dif)>0.0))
  176.     error, "***WARNING*** SVdec returned increasing singular values";
  177.   achk= u(,+) * (s*v)(+,);
  178.   sabs= max(abs(s));
  179.   err= max(abs(a-achk))/sabs;
  180.   if (err>1.e-9) {
  181.     write, "***WARNING***  SVdec decomposition doesn't check";
  182.     write, "   max relative error is "+pr1(err);
  183.   }
  184.  
  185.   s= SVdec(a, u, v, full=1);
  186.   if (anyof(s(dif)>0.0))
  187.     error, "***WARNING*** SVdec returned increasing singular values";
  188.   uu= u(,1:min(m,n));
  189.   vv= v(1:min(m,n),);
  190.   achk= uu(,+) * (s*vv)(+,);
  191.   sabs= max(abs(s));
  192.   err= max(abs(a-achk))/sabs;
  193.   if (err>1.e-9) {
  194.     write, "***WARNING***  SVdec decomposition doesn't check";
  195.     write, "   max relative error is "+pr1(err);
  196.   }
  197.   err= max(abs(u(,+)*u(,+)-unit(m)))+max(abs(v(,+)*v(,+)-unit(n)));
  198.   if (err>1.e-9) {
  199.     write, "***WARNING***  SVdec decomposition not orthogonal";
  200.     write, "   max relative error is "+pr1(err);
  201.   }
  202. }
  203.  
  204. func testSV(n)
  205. {
  206.   a= random(n,n);
  207.   b= random(n);
  208.   SVcheck,a,b,SVsolve(a,b), "1D";
  209.   b2= random(n);
  210.   x= SVsolve(a,[b,b2])
  211.   SVcheck,a,b, x(,1), "2D(1)";
  212.   SVcheck,a,b2, x(,2), "2D(2)";
  213.   x= SVsolve(a,transpose([b,b2]), which=2)
  214.   SVcheck,a,b, x(1,), "2D(1)/which";
  215.   SVcheck,a,b2, x(2,), "2D(2)/which";
  216. }
  217.